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A Geometric  Approach  for  Knot  Selection  in 
Convexity-Preserving  Spline  Approximation 

R.  Morandi,  D.  Scaramelli,  and  A.  Sestini 


Abstract.  A geometric  approach  is  proposed  for  selecting  the  knots  used 
in  a parametric  convexity-preserving  B-spline  approximation  scheme.  The 
approach  automatically  gives  the  necessary  information  about  the  shape 
suggested  by  the  data  which  may  be  exact  or  not. 


§1.  Introduction 

In  many  different  fields  such  as  medicine,  physics,  engineering  and  computer 
graphics  the  amount  of  data  obtained  through  experimental  and/or  statistical 
surveys  is  very  large.  Consequently,  the  selection  of  a suitable  small  set  of 
knots  becomes  an  indispensable  step  within  any  efficient  spline  approximation 
scheme. 

In  this  paper,  large  sets  of  exact  and  non-exact  data  points  are  approxi- 
mated by  means  of  a spline  approximation  scheme.  So  a knot  selection  strat- 
egy is  necessary  and  this  is  provided  by  a geometric  approach.  In  detail,  the 
proposed  approach  is  based  on  some  weights  suitably  associated  to  the  data 
points  and  directly  computed  from  them.  In  addition,  the  method  automati- 
cally defines  the  shape  suggested  by  the  data,  here  assumed  planar  and  exact 
or  not.  Furthermore,  the  shape  constraints  for  the  approximating  curve  can 
be  obtained  in  order  to  reproduce  the  desired  behaviour. 

Several  other  approaches  have  been  studied,  such  as  the  knot  removal 
methods  [5,6],  to  reduce  the  number  of  parameters  involved  in  an  approxima- 
tion problem.  Interesting  results  are  already  available  even  for  constrained 
approximation  [1],  However,  the  approach  introduced  here  differs  from  those 
because  it  does  not  reduce  an  initial  large  set  of  knots,  but  directly  computes  a 
suitable  set.  Furthermore,  it  does  not  require  the  starting  approximation  with 
many  knots  used  in  [1,5,6]  for  the  weight  definition  (solving  a minimization 
problem  for  each  weight). 

The  proposed  strategy  has  been  tested  on  several  examples  for  open  and 
closed  curves,  and  for  exact  and  non-exact  data  points.  The  approximating 
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curve  is  obtained  by  means  of  a convexity-preserving  B-spline  approximation 
scheme.  The  goodness  of  fit  of  the  approximation  has  been  estimated  mea- 
suring the  mean  square  distance  of  the  data  points  from  the  resulting  curve. 

The  outline  of  the  paper  is  as  follows.  The  problem  and  the  method 
are  presented  in  the  next  section.  The  strategies  used  to  define  the  shape 
suggested  by  the  data  and  to  select  the  knots  are  introduced  in  Section  3, 
and  they  are  given  in  more  detail  in  the  Appendix.  Finally,  in  Section  4 the 
numerical  results  obtained  for  four  sets  of  data  points  are  given. 

§2.  The  Problem  and  the  Method 

Let  Pi  £ 1R2,  i = 0, . . . , n,  be  exact  or  non-exact  data  points,  with  n very  large, 
and  let  T = {t*  £ 1R,  i = 0, . . . , n}  be  an  assigned  set  of  associated  strictly 
increasing  parameter  values.  For  non-exact  data,  let  d be  a positive  assigned 
quantity  so  that  ||Pj  — P\ ||2  < d,i  = 0,...,n,  where  Pf  is  the  (unknown) 
exact  data  point  corresponding  to  P,.  We  observe  that,  for  simplicity’s  sake, 
the  same  maximum  error  value  is  assumed  for  all  the  data  points. 

It  is  well  known  that,  if  the  data  sets  are  very  large  and  in  particular  the 
data  are  non-exact,  the  use  of  an  approximation  scheme  is  the  only  reasonable 
approach  to  construct  a curve  with  the  desired  behaviour.  Thus,  here  the  aim 
is  to  first  define  the  convexity  constraints  suggested  by  the  data,  and  then  to 
give  a strategy  for  selecting  a suitable  small  number  of  knots  to  construct  a 
convexity-preserving  least-square  B-spline  approximating  curve. 

Let  Njk(t),j  — 1 — k, . . . ,nr  — 1,  be  the  usual  B-splines  of  order  k [4]  de- 
fined with  an  extended  knot  vector  0*  = {ri_k, . . . , To, . . . , Tnr, . . . , Tnr+k-i}- 
Thus,  we  can  introduce  the  B-spline  representation  of  a spline  curve 

nr— 1 

C(t)  = QjNJk{t),  t £ [to,  Tnr],  (1) 

j~l  — k 

where  Qj,j  = 1 — k, . . . ,nr  — 1 are  de  Boor  control  points. 

The  problem  can  be  divided  into  three  sub-problems.  The  definition  of 
the  shape  suggested  by  the  data  (that  is  the  determination  of  the  convex- 
ity changes  required  to  the  approximating  curve),  the  selection  of  the  knots, 
and  finally  the  construction  of  the  convexity-preserving  least-square  B-spline 
approximating  curve. 

In  particular,  the  shape  suggested  by  the  data  is  obtained  through  the 
procedure  called  “Shape  Determination”  (SD)  described  in  the  Appendix.  SD 
uses  some  coefficients  = 0,  ...,n,  suitably  associated  with  the  data  to 
establish  in  which  parameter  values  zero  curvature  is  required,  and  to  deter- 
mine the  curvature  sign  in  the  interval  [To,rnr],  As  the  planar  case  is  here 

considered,  the  curvature  is  defined  as  the  function  pit)  — C'{t)xC(t)  wjiere 

V ' 110(0111 

v x z = v\  ■ Z2  — i»2  • 2i,  Vv,z  £ 1R2).  Thus,  the  procedure  “Knots  Selection” 
(KS)  described  in  the  Appendix  selects  the  knot  vector  0 = {to,  . . . , rnr } by 
using  the  weights  ra,  = ju,  |,  i = 0,  0*  is  defined  as  the  corresponding 
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extended  knot  vector,  taking  into  account  whether  the  approximating  curve 
is  open  or  closed. 

Finally,  the  last  step  is  realized  through  the  solution  of  a constrained  para- 
metric least-squares  problem  as  a general  constrained  optimization  problem. 
In  fact,  as  the  general  parametric  case  is  considered,  the  objective  function 
E?=o  II pi  - C'(t»)ll 2 = E1U  ll-P*  - EjET-fc  QjNik(ti)\\l  is  quadratic  in  the 
unknowns  Qj,j  = 1 — 1,  but  the  convexity  constraints  are  nonlinear. 


§3.  Data  Shape  Determination  and  Knot  Selection  Strategy 

The  SD  procedure  defines  the  shape  suggested  by  the  data,  and  KS  selects 
the  knots  for  constructing  the  B-spline  curve.  They  are  presented  in  the  Ap- 
pendix, but  are  commented  upon  here.  Concerning  the  shape  determination, 
SD  computes  the  set  U = {«o, . . . , un}  whose  sign  variations  are  the  cur- 
vature sign  variations  required  for  the  approximating  curve.  If  exact  data 
are  considered,  |it,|  is  the  reciprocal  of  the  radius  of  the  circle  which  passes 
through  the  points  P;_i,  Pi,  Pj+i,  and  its  sign  is  that  of  A*  = \det(Li,  Li+i), 
where  Li  = Pi  — Pf_i.  On  the  other  hand,  the  sign  of  A*  is  a reliable  geo- 
metric information  in  the  case  of  non-exact  data  only  if  the  condition  (2)  of 
the  theorem  given  in  this  section  holds.  Thus,  the  definition  of  U is  suitably 
modified  for  non-exact  data,  using  an  input  tolerance  told  and  considering  the 
result  given  in  the  theorem  below.  In  this  case  Ui  is  defined  using  the  circle 
through  the  points  Pi, , Pi  and  PTi , where  P/t  and  PTi  are  suitably  selected 
points.  In  detail,  < i - 1,  r<  > * + 1,  Efc=(,  ||Pfc+i  ~ Pfclh  < told  • Lp  and 
EfcEi1  ||Pfc+i  — Pfclh  < told  ■ Lp,  where  Lp  is  the  length  of  the  polygonal 
joining  the  data  points.  SD  computes  the  set  0s  = {tq,  . . . , t*s  } c T and 
the  set  E = [oq,  . . . ,<rns-i},  where  &s  is  such  that  Tq  = to,  r®s  = tn  and 
zero  curvature  is  required  at  each  r*,  i = 1, . . . , ns  — 1.  The  desired  curvature 
sign  between  r*  and  r®+1  is  given  by  cq  equal  to  -1  or  0 or  1.  We  observe 
that,  for  d ^ 0,  in  the  procedure  w,  ^ 0 is  assumed  to  imply  the  existence  of 
ki  £ {i  — 2 ,i  — 1,  z}  such  that  ■ Ui  > 0,  ■ u*  > 0 and  u^+2  • tq  >0. 

This  hypothesis  seems  to  be  quite  reasonable  as  n is  considered  very  large. 

Theorem  1.  Let  Pi  E 3R2  for  i = 0, . . . , n be  assigned  non-exact  data  points, 
and  let  d be  a small  positive  assigned  quantity  such  that  1 1 Pi  — P®  1 12  < d,  i = 
0, . . . ,n,  where  P\  is  the  (unknown)  exact  data  point  corresponding  to  Pi.  If 
d < \ mini=ii...in  ||Lj||2  with  L,  ~ Pi  — P;_i  and  the  condition 


, 11, 

\\Li\\2  + ||Li+1||2  > 4a- 


i = 1, . . . ,n  — 1, 


(2) 


holds,  then 

Nei  -Ni  >0,  * = l,...,n  — 1, 

where  Ni  = L,  A Lj+i,  N ® = L®  A Lei+1,  Le{  = P|  — P\_  j and  the  symbols 
” A ” and  denote  the  usual  vector  and  scalar  product,  respectively. 

Proof:  We  can  write  P\  = P,  + with  0 < £i  < d and  ||tq||2  = 1.  Then 
we  have  Lf  = i,  +£jWi  -ej-iw,-!  and  IV * = + zit  where  zt  = ZiVif\Li+i  - 
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A Li+ 1 + £i  + l Li  A t>i  + i + £j£j-f  1 Vi  A Vi+1  — £j_i£i  + iVi_l  A Vi  + 1 — 
EiLi  A V;  +£j-i£iVj_i  A V;.  Thus,  we  can  write  TV-  = (1  + /(£i_i,£j, £i+i))Wj, 
where  /(£j_i,£j,£j+i)  = As  a consequence,  the  assertion  holds  if 

(1  + /(£j_i,£i,ej+i))||iVj||^  > 0,  that  is  if  /(£i_i,£j,£i+i)  > -1.  Now  with 
some  algebra  the  following  inequality  can  be  easily  obtained: 

/(e«-i,ei,e«+i)  > - ||jy.||2 Pdl-^lh  + ll-^i+ilb)  + 3d].  (3) 


So,  as  d < ^ mini=i  ..  n | |Z/z- 1 12 , d < ll-^lb-HlT.-i-ilb  . fjjen,  the  inequality  (3) 
and  the  hypothesis  (2)  imply  that  /(£,_!,£;, £i+i)  > - ^ (xOl^lk  + 
H2))  > — 1)  thus  proving  the  theorem.  □ 

Concerning  the  knot  selection,  in  the  KS  procedure  the  knot  vector  OcT 
is  initialized  with  0.5-  Then,  for  each  interval  [r?_i, rj’], j = l,...,ns,  it  is 
checked  if  other  parameters  must  be  inserted  in  the  knot  vector  using  weights 
Wi  — |n;|,  i = 0, . . • , n.  More  precisely,  a parameter  value  U &T  D (rj_1,  T-)  is 
inserted  in  0 if  one  of  the  following  conditions  holds:  either  the  corresponding 
weight  wi  is  big  enough  and  Pi  is  far  enough  from  all  the  other  data  related  to 
the  parameters  previously  introduced  between  r?_  j and  r®,  or  the  correspond- 
ing weight  wi  is  not  big  enough  but  Pi  is  too  far  from  them.  For  choosing 
reasonable  values  for  the  tolerances  used  in  the  previous  consideration,  de- 
noted as  tolw,tol,n  and  told 2,  it  is  assumed  that  the  distances  are  relative  to 
an  approximated  curve  length  and  the  weights  are  relative  to  the  maximum 
weight.  The  parameter  values  between  rj '_i  and  r®  are  ordered  according  to 
a decreasing  order  of  the  corresponding  weights. 

§4.  Numerical  Results 

Four  numerical  tests  are  presented  to  analyze  the  performance  of  the  approach. 
For  all  the  considered  tests  the  parameter  values  t*, i = 0, . . . ,n,  have  been 
computed  with  the  chord-length  approach  and  a scaling  such  that  0 = <0  < 
t !<•••<  tn_i  < tn  = 1.  The  approximated  curve  length  L required  as  input 
by  the  KS  procedure  has  been  computed  as  the  length  of  the  piecewise  linear 
interpolant  of  all  the  data  points  if  d = 0 and  of  a suitably  selected  subset  of 
them  if  d > 0. 

A sequential  quadratic  programming  method  [2]  is  used  to  construct  the 
approximating  curve  by  means  of  the  routine  CONSTR  of  the  Optimization 
toolbox  of  the  Matlab  package  [3].  The  set  of  control  points  used  to  start  is 
chosen  as  the  set  of  data  points  corresponding  to  the  selected  knots.  Con- 
cerning the  constraints,  as  k = 4 has  been  used  in  the  experiments,  four 
consecutive  control  points  are  required  to  generate  a convex  polygon  if  con- 
vexity is  looked  for  in  the  corresponding  curve  segment.  In  addition  three 
sequential  collinear  control  points  are  required  if  zero  curvature  is  asked  at 
the  corresponding  knot. 
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For  each  test,  a figure  is  given  showing  the  corresponding  approximating 
curve  on  the  left,  and  the  related  curvature  plot  on  the  right.  In  all  the 
figures  the  data  points  Pq,  ...,  Pn  are  denoted  with  the  symbol  “ ■ the  points 
associated  with  the  knots  with  the  symbol  “ o ” , and  the  points  corresponding 
to  the  knots  belonging  to  ©5  with  the  symbol  “ ® 

The  results  are  summarized  in  Table  1,  where  n + 1 is  the  number  of  data 
points,  told  is  the  input  tolerance  used  in  the  procedure  SD,  tolw,toldi  and 
told2  are  the  input  tolerances  used  in  the  procedure  KS  and  nr  + 1 is  the 
cardinality  of  the  knot  vector  0.  The  symbol  % denotes  the  percentage  ratio 
(nr  + 1 )/(n  + 1).  To  estimate  the  goodness  of  fit  of  the  approximation,  the 
Mean  Square  Distance  (MSD)  of  the  data  points  P*,i  = 0,  from  the 

approximating  curve  is  also  given  in  the  table. 


Tab.  1.  Results  of  the  tests. 

In  the  first  test,  data  are  considered  exact  ( d = 0),  while  in  the  other 
tests  they  are  non-exact. 

Test  1 relates  to  a set  of  285  exact  data  points  that  represent  the  alpha- 
bet capital  letter  “ D ”.  In  this  case  only  a curvature  sign  variation  to  the 
approximating  curve  is  required,  as  the  curvature  plot  shows. 

In  Test  2,  257  non-exact  data  points  are  considered.  They  have  been 
obtained  by  introducing  a simulated  random  perturbation  with  d = 0.4  on  the 
ordinates  of  the  points  (xf,yf),i  = 0, . . . , 256  defined  as  x\  — -8  + dx-i,  yf  = 
1251^f'^  = 0,  ...,256,  where  R,  = sqrt(  2(xf)2)  + eps,  with  eps  denoting 
the  round-off  error  and  dx  = 0.0625.  We  can  observe  that  the  simulated  error 
does  not  preserve  the  symmetry  of  the  data,  and  therefore  the  selected  knots 
are  not  symmetric. 


0 0.2  0.4  0.6  0.8  1 


2.  Test  2:  On  the  left  the  approximating  curve,  and  on  the  right  the  related 
curvature  plot. 


In  Test  3 the  error  in  the  126  data  is  obtained  by  introducing  a simulated 
random  perturbation  with  d = 0.12  on  both  the  coordinates  of  the  points 
i = 0, . . . , 125  defined  as  sf  = sin{A-n  ■ dt  ■ i),  yf  = cos(2-7r  • dt  ■ i),  i = 
0, . . . , 125,  where  dt  = 0.008.  In  this  case,  among  the  data  there  are  sequences 
of  quasi-collinear  points,  and  we  can  observe  that  the  approximating  curve 
has  corresponding  almost  straight  line  segments. 

An  application  to  an  engineering  problem  is  presented  in  Test  4.  The  244 
data  points  are  derived  from  measurements  effected  in  the  Power  Station  lo- 
cated in  Seraing  (Belgium).  The  measurements  are  related  to  the  active  power 
of  the  alternator  in  the  Central,  observed  on  03/01/1997  between  6:30:00  and 
7:00:00.  In  this  case  the  maximum  value  of  the  error  is  d = 1. 

It  should  be  noted  that,  the  user  needs  to  work  in  an  interactive  way  for 
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Fig.  3.  Test  3:  On  the  left  the  approximating  curve,  and  on  the  right  the  related 
curvature  plot. 


Fig.  4.  Test  4:  On  the  left  the  approximating  curve,  and  on  the  right  the  related 
curvature  plot. 

selecting  suitable  values  of  the  input  tolerances.  Obviously,  a previous  check 
of  the  data  shape  and  of  the  data  distribution  is  of  great  help. 

§5.  Appendix 


SD  Procedure  (Shape  Determination) 


Input:  V = {P0,  . . . , P„},  d,  told , T = {t0, . . . , tn} 

• Define  two  auxiliary  suitable  data  points  P_j  and  Pn+i 
if  d = 0 

for  i = 0, . . . , n 

ll-MHI-fo+ilIHlV'ilh 


•Ui  = 
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with  Li  = Pi  - Pi- 1,  Vi  - Pi+ 1 - Pi- 1,  A,  = \det(Li,Li+i) 

end 

else 

•Lp  = ElZl\\Pk+i-Pk\\2 

for  i = 0, . . . , n 

•levTi  = maxJ=i+li...in+i{j  : ||P*+i  - f’ifclh  < told  ■ Lp } 

•levh  = {j  : ll-Pfc+i  - -Pfclk  < told  ■ Lp} 

end 

for  i = 0, . . . , n 

• determine  the  biggest  U 6 {— 1, . . . , i — 1}  and  the 
smallest  r*  € {i  + 1, . . . , n + 1}  such  that  (2)  holds 
replacing  Pj+i  with  Pr<  and  P;_i  with  P i{ 

if  li  < levi{  or  r i > levri 

•Mj  = 0 

else 

#M.  = JAi . 

1 IlI-db'llii-Hlh'llVjIb 

where  Li  =Pi~  Pi,,  Li+1  = Pr>  - Pi,  Vi  = PTi  - Pu, 

Ai  = r^d^t^Li,  Li-i-i} 

end 

end 

end 
•i  = 1 
•Os  = {to} 
while  i < n — 1 
if  Mi  = 0 

•Os  = Os  U {<i-i} 

• determine  0<id<n  — i — 1 such  that 

Mi  = . . . = Mi+ij  = 0 and 
Mi+id+i  ^ 0 or  id  = n — i — 1 
•Os  = ©s  U {t|^2i±i!j } 

•Os  = ©s  U {<i+id+i} 

•i  = i + id  + 2 
elseif  Mi  • Mi+i  < 0 

if  |Mi|  < |Mi+1| 

•Os  = ©s  U (tj 
•i  = i + 1 

else 

•Os  = ©s  U {<i+i} 

•i  = i + 2 

end 

else 

•i  = i + 1 

end 

end 

•Os  = ©s  U {t„}  d=  {r0s, . . . ,t*s} 
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for  j = 0, . . . , ns 

•s(j)  =i  where  r?  = U 

end 

for  j = 0, . . . ,ns  — 1 

if  s(j  + 1)  = s(j)  + 1 

•Oj  = 0 

else 

•Oj  = 

end 

end 

Output:  U = {u0,...,un},  Qs  = {t0s,...,t*s},  E = {cr0,  • • • , o^-i}- 


KS  Procedure  (Knot  Selection) 

Input: 

V = {P0,  • • • , Pn},  T = {to,  ■ ■ ■ , t«}, 

U = {u0,  ©5  = {Tq  ) • • • i Tns}> 

L i tolw , t() l , told 2 {toldl  ^ told2 ) 

•Wi  = f = 0, . ..,n 

•0  = ©s 

•»mai  = max{«;o, . . . , w„} 

for  j = 0, . . . , 7is 

•s(j)  = i where  t?  = U 

end 

for  j = 1, . . . ,ns 

•Sj  = s(j)  - s(j  - 1)  + 1 

• let  {*i, } be  the  index  permutation  of  {s(j  — 1), , s(j)} 
such  that  the  weights  , . . . , are  in  decreasing  order 

*®J  = {^s(j-l)i  ^s(;')} 

= {-PsO-l),^)} 
for  fc  = 1, . . . , Sj 

if  ^>to1”  and  v^€Pr3.  ^~Pr"2  >toldl 
•©J  = © j © {tik  } 

•VRj  =VRj  U{PiJ 

elseif  < tolw  and  VPr  € VRl  l|J*it=PilU  > told2 

UJTnax  J Lj 

•Qj  = 0jU{tj,.} 

•PRj  = 'PRj  U{-PiJ 

end 

end 

•©  = 0 U Qj 

end 

Output:  0 = {r0 ,...,Tnr} 
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